Imagine you are a newbie surveyor sent to the field. You only have two instruments, your brain and your phone. Your boss, (and he's harsh) requires that you create a terrain profile of your journey from Point A to Point B. And the deadline is just 7 days away. What do you do? Most of today's smartphones can download a GPS app, which also subsequently calculate the altitude of an area. Problem one sorted. You collect the GPS locations with their corresponding elevation readings and call it a day. Problem two: how to plot the terrain profile in python. You remember that Qgis has a very simple tool called the Profile tool, but then a flashing thought reminds you that you were teaching yourself how to code :(. Practice makes perfect, and you have to perfect the art of your coding practice.
As usual, you go to your laptop and as always, you begin by loading the usual tools.
# Import the needed packages for data science and dataframes
import pandas as pd
import geopandas as gpd
import os
Things are proceeding well, you have already loaded the necessary packages. Next, you want to call in the csv file containing your GPS readings. That's also an easy task. You actually do it as you engage in an impromptu phone call.
# Import the CSV containing the GPS coordinates and their altitude
path = "E:/documents/gis800_articles/jupyter/terrain_profile"
gps = pd.read_csv(os.path.join(path, 'gps.csv'))
gps.head()
| Point | Latitude | Longitude | Accuracy(m) | Altitude(m) | Description | |
|---|---|---|---|---|---|---|
| 0 | 1 | -1.390684 | 36.769671 | 4 | 1709 | Junction at Magadi road |
| 1 | 2 | -1.398675 | 36.790664 | 9 | 1682 | Start of murram road |
| 2 | 3 | -1.397448 | 36.798129 | 3 | 1660 | Bend before reaching the bridge crossing river... |
| 3 | 4 | -1.397348 | 36.798175 | 5 | 1655 | a few metres to the bridge |
| 4 | 5 | -1.399695 | 36.799621 | 3 | 1674 | A group of zebras were sighted grazing in an o... |
The above steps were a no brainer. You even commend yourself for being a good programming learner. Small wins lead to big wins. You actually look forward to a slightly harder task of converting the above csv table to a dataframe. This is because as per your programming lessons so far, a geodataframe has to be created from a dataframe if you want to work with tabular data for geospatial analyses. This is another exciting exercise, and because of its simplicity, you do it in a whizz.
# Convert to a dataframe
gps_df = pd.DataFrame(gps)
gps_df.head()
| Point | Latitude | Longitude | Accuracy(m) | Altitude(m) | Description | |
|---|---|---|---|---|---|---|
| 0 | 1 | -1.390684 | 36.769671 | 4 | 1709 | Junction at Magadi road |
| 1 | 2 | -1.398675 | 36.790664 | 9 | 1682 | Start of murram road |
| 2 | 3 | -1.397448 | 36.798129 | 3 | 1660 | Bend before reaching the bridge crossing river... |
| 3 | 4 | -1.397348 | 36.798175 | 5 | 1655 | a few metres to the bridge |
| 4 | 5 | -1.399695 | 36.799621 | 3 | 1674 | A group of zebras were sighted grazing in an o... |
You love the progress. All the above done in less than 10 minutes. You just want to confirm that what you've got is indeed a dataframe. The feeling of transforming a csv to a dataframe should be enough reward for the coding labour so far.
# Check the class of gps_df
type(gps_df)
pandas.core.frame.DataFrame
Now, the following step of converting the dataframe to a geodataframe is not too intuitive. It requires more parameters and arguments than the above four coding steps. However, based on your (sometimes painful) programming tutorials, you have been taught the art of searching for solutions from the web. You embark on this, and you learn that geopandas has the GeoDataFrame method for precisely this function. It's another walkover.
However, there is a twist. The GeoDataFrame method requires we define the CRS. So you quickly search for which package does this, and then call it into Jupyter, your working platform.
# Import the package for CRS manipulations
from pyproj import CRS
# Now convert to a geodataframe
gps_points = gpd.GeoDataFrame(gps_df,
geometry=gpd.points_from_xy(gps_df['Longitude'], gps_df['Latitude']),
crs=CRS.from_authority('EPSG', '4326')) # Assigns the CRS EPSG: 4326 to the geodataframe
gps_points.head()
| Point | Latitude | Longitude | Accuracy(m) | Altitude(m) | Description | geometry | |
|---|---|---|---|---|---|---|---|
| 0 | 1 | -1.390684 | 36.769671 | 4 | 1709 | Junction at Magadi road | POINT (36.76967 -1.39068) |
| 1 | 2 | -1.398675 | 36.790664 | 9 | 1682 | Start of murram road | POINT (36.79066 -1.39868) |
| 2 | 3 | -1.397448 | 36.798129 | 3 | 1660 | Bend before reaching the bridge crossing river... | POINT (36.79813 -1.39745) |
| 3 | 4 | -1.397348 | 36.798175 | 5 | 1655 | a few metres to the bridge | POINT (36.79817 -1.39735) |
| 4 | 5 | -1.399695 | 36.799621 | 3 | 1674 | A group of zebras were sighted grazing in an o... | POINT (36.79962 -1.39969) |
You clap yourself for the progress so far. Once again, you want to exercise your fingertips once more by confirming the above table is indeed a geodataframe.
# Check the CRS of our geodataframe
gps_points.crs
<Geographic 2D CRS: EPSG:4326> Name: WGS 84 Axis Info [ellipsoidal]: - Lat[north]: Geodetic latitude (degree) - Lon[east]: Geodetic longitude (degree) Area of Use: - name: World. - bounds: (-180.0, -90.0, 180.0, 90.0) Datum: World Geodetic System 1984 ensemble - Ellipsoid: WGS 84 - Prime Meridian: Greenwich
Okay, so far so good. You were dying to draw a terrain profile. Wasn't that the prime aim in the first place? Too much aimless coding will not replace that your boss wanted a terrain profile of your fieldwork in less than 7 days (minus a few hours by now). You search online. From the videos on YouTube such as this one here that last 45 minutes, and as if that is not enough, they have a sequel to complete the job, drawing a simple profile is now seems scary. "Why should something so simple take so long?", you ask yourself. You have the temptation of using Qgis, but it will be like a blow to your coding ego. To preserve your machoism, you remember that Plotly is a package for creating cool graphs. Why not use that to create an elevation profile? Again, you first need to rope in the relevant packages.
# Plotly is a tool for drawing interactive graphs
import plotly.express as px
import kaleido # Useful in saving plotly graphs to a directory
# Plot the terrain profile
fig = px.line(data_frame=gps_points, x='Point', y='Altitude(m)', text='Description')
fig.update_traces(textposition='top center')
fig.update_layout(
height=800,
title_text='Terrain Profile from Ongata Rongai to Tasia, Kajiado County'
)
fig.show()
fig.write_image(os.path.join(path, 'terrain.jpg'), width=1500, height=900)
That's it! Now to just confirm that indeed these coordinates fall on the right places in the real world, you have to plot them overlying a basemap. A basemap with your GPS points will not only confirm that the datum and their positioning is correct, but it may also iron out any trust issues your supervisor may have, such as suspecting you did this via Google map instead of going to the field. Did someone mention basemaps? As per your research, contextily is a wonderful package for wonderful maps.
# Import contextily
import contextily as cx
import matplotlib.pyplot as plt
Now that you have you have loaded the requisite packages, time for the ultimate delivery: a basemap with your GPS coordinates and their corresponding elevations labelled.
fig, ax = plt.subplots(figsize=(15, 9))
gps_points_map = gps_points.plot(ax=ax)
for alt, lon, lat in zip(gps_points['Altitude(m)'], gps_points['Longitude'], gps_points['Latitude']):
gps_points_map.annotate(alt, xy=(lon, lat), xytext=(-4, 2), textcoords='offset points', color='white')
cx.add_basemap(ax, crs=CRS.from_authority('EPSG', '4326'), source=cx.providers.Esri.WorldImagery)
Let's suppose you would also like to plot the path you followed as you walked with your phone whilst collecting the GPS coordinates. This also sounds easy, but is far from it. Congrats for the job so far, but we take over from here. We shall begin by converting the geodataframe to a LineString. To work with linestrings, we use the shapely.geometry package.
# Import the package to work with linestrings
from shapely.geometry import LineString
Let's create afor loop that converts a geodataframe to a linestring.
# Extract list of coordinates to create linestring
geom_list = []
for xy in zip(gps_points['Longitude'], gps_points['Latitude']):
geom_list.append(xy)
print(geom_list)
[(36.76967061, -1.39068412), (36.79066426, -1.3986755), (36.79812947, -1.3974481), (36.79817471, -1.39734767), (36.79962135, -1.39969483), (36.80488978, -1.40220242), (36.81747036, -1.40132149), (36.8218044, -1.40157005), (36.82966621, -1.40291262), (36.8349445, -1.40457771), (36.82987579, -1.40756392), (36.82985959, -1.40834375)]
geom_path = LineString(geom_list)
geom_path
Alright. We have a LineString. Check. But can it be plotted on a basemap? A resounding 'No'. As always, we went the extra mile of checking it out ourselves so save your energy for something else. Apparently, it appears that LineStrings are only good for displaying the geometry of an object over a Cartesian plane, not the real earth. To drive the point home, the LineString has no geographical reference whatsoever. What you see is just a figure following what may read as a list of coordinates, but in actual sense they are being read as cartesian coordinates. If the LineString had a geospatial reference, we could plot it over a contextily basemap straight away and end it here but this is far from the truth! Because of this shortcoming, we need a more portable format. Good ol' shapefile(s) never disappoint. If you have a genius mind (of course you do) you can already perceive all that is needed is to convert the LineString to a shapefile. Yes, you're right. However, it is not so straightforward. Luckily, some clever fellows at Hatari Labs have already charted out a way for us to do this kind of thing, thus saving us from a near mental breakdown.
As per their workflow, they transform a geodataframe--in our case gps_points--to a shapefile. Therefore let's go back to our roots. Call gps_points back to duty.
# Remind yourself of the geodataframe structure
gps_points.head()
| Point | Latitude | Longitude | Accuracy(m) | Altitude(m) | Description | geometry | |
|---|---|---|---|---|---|---|---|
| 0 | 1 | -1.390684 | 36.769671 | 4 | 1709 | Junction at Magadi road | POINT (36.76967 -1.39068) |
| 1 | 2 | -1.398675 | 36.790664 | 9 | 1682 | Start of murram road | POINT (36.79066 -1.39868) |
| 2 | 3 | -1.397448 | 36.798129 | 3 | 1660 | Bend before reaching the bridge crossing river... | POINT (36.79813 -1.39745) |
| 3 | 4 | -1.397348 | 36.798175 | 5 | 1655 | a few metres to the bridge | POINT (36.79817 -1.39735) |
| 4 | 5 | -1.399695 | 36.799621 | 3 | 1674 | A group of zebras were sighted grazing in an o... | POINT (36.79962 -1.39969) |
We will create the format in which our attributes shall be stored. For now, we shall give it the buzzword schema. The dictionary format that goes into the schema will actually be the blueprint for how our shapefile's attribute table shall appear. Just consider the schema as a framework.
# define schema
schema = {
'geometry':'LineString',
'properties':[('Name', 'str')]
}
We will also import the fiona package to help us in creating a shapefile.
# Import fiona
import fiona
Let's create our shapefile. Note that the schema parameter will arrange the attributes as per the dict format we specified.
lineShp = fiona.open(os.path.join(path, 'shapefile', 'path.shp'), mode='w', driver='ESRI Shapefile',
schema = schema, crs = "EPSG:4326")
To create a shapefile that will resemble our LineString, we have to extract the coordinates from geom_path in the same order, that is, as the crow flies. The for loop for this job shall iterate through every row and collect the coordinates to be imputed to our shapefile --lineShp. Remember the lineShp has no attributes in it whatsoever, but we had designed a framework (called schema) on how they should appear. The for loop will just extract these and fit them into the schema template.
#get list of points
xyList = []
rowName = ''
for index, row in gps_points.iterrows():
xyList.append((row['Longitude'], row['Latitude']))
rowName = 'shapefile' # assign a string value 'shapefile' to each row
print(xyList)
print(rowName)
[(36.76967061, -1.39068412), (36.79066426, -1.3986755), (36.79812947, -1.3974481), (36.79817471, -1.39734767), (36.79962135, -1.39969483), (36.80488978, -1.40220242), (36.81747036, -1.40132149), (36.8218044, -1.40157005), (36.82966621, -1.40291262), (36.8349445, -1.40457771), (36.82987579, -1.40756392), (36.82985959, -1.40834375)] shapefile
All right. Time to write our attributes to the shapefile. We will create a dictionary called lineDict that follows the schema format we had defined earlier. We will define the geometry format of our shapefile as LineString and for the coordinates, they shall come from the xyList list. For the properties, only a word is needed. It sounds confusing, but remember this is a line vector we are creating. It only has one row(!) unlike a point shapefile which can have hundreds. Thus, only one string value is needed to get into the Name attribute property. We chose shapefile. This insertion of a Name property may sound trivial, but omit it and it's impossible to create the shapefile. Just try.
#save record and close shapefile
lineDict = {
'geometry' : {'type':'LineString',
'coordinates': xyList},
'properties': {'Name' : rowName},
}
lineShp.write(lineDict)
#close fiona object
lineShp.close()
Now let's load our shapefile and plot it on a basemap without further ado.
# Read shapefile
shapefile = gpd.read_file(os.path.join(path, 'shapefile', 'path.shp'))
shapefile.plot()
<AxesSubplot: >
Lastly, we plot the path that we followed as we were conducting our fieldwork. Truth of the matter is, we wished we could plot the altitude points and our path on the same contextily map, but searching for a way to do this was fruitless. Here is your path anyway.
# Plot your path
fig, ax = plt.subplots(figsize=(15, 9))
shapefile_map = shapefile.plot(ax=ax)
cx.add_basemap(ax=ax, crs=CRS.from_authority('EPSG', '4326'), source=cx.providers.Esri.WorldImagery)
That's it!
Just as one can plot a terrain profile in Qgis using the Profile tool, the same can be accomplished in Python with just altitude values present. Where the elevation values are missing, a Digital Elevation Model (DEM) will have to be present so that the height values at the Lon-Lat coordinates are extracted and imputed to the table with coordinates. In this article, we have explained the process of creating a terrain profile from a csv file. The csv file is first transformed to a geodataframe, and by use of the plotly package, a terrain profile can be deduced. In a possible future exercise, we shall investigate creation of a terrain profile using a DEM.